Simulating Stochastic Inertial Manifolds 
by a Backward- Forward Approach * 



Xingye Kan^'^, Jinqiao Duan^'^, loannis G. Kevrekidis^ and Anthony J. Roberts^ 

1. Department of Applied Mathematics 
Ilhnois Institute of Technology 
Chicago, IL 60616, USA 
E-mail: xkan@hawk.iit.edu, duan(Qliit.edu 

2. Institute for Pure and Applied Mathematics, UCLA 
Los Angeles, CA 90095, USA 

3. Program in Applied and Computational Mathematics 

& 

Department of Chemical and Biological Engineering 
Princeton University 
Princeton, NJ 08544, USA 
E-mail: yannis@princeton.edu 

4. School of Mathematical Sciences 
University of Adelaide 
Adelaide, Australia 
E-mail: anthony.roberts@adelaide.edu.au 

June 18, 2012 



Abstract 

A numerical approach for the approximation of inertial manifolds of 
stochastic evolutionary equations with multiplicative noise is presented 
and illustrated. After splitting the stochastic evolutionary equations 
into a backward and a forward part, a numerical scheme is devised 
for solving this backward-forward stochastic system, and an ensemble 
of graphs representing the inertial manifold is consequently obtained. 
This numerical approach is tested in two illustrative examples: one is 
for a system of stochastic ordinary differential equations and the other 
is for a stochastic partial differential equation. 



*This work was partly supported by the NSF Grants 1025422, 0923111 and 0731201, 
by the US DOE grant DESC0002097, and by the NSFC grants 10971225 and 11028102. 



1 



Key words Invariant manifolds, inertial manifolds, random dy- 
namical systems, backward and forward stochastic differential equa- 
tions, stochastic partial differential equations, numerical schemes 



1 Introduction 

The concept of inertial manifolds for deterministic partial differential equa- 
tions was introduced in 1980s [51 [TB]. These manifolds are finite dimensional 
invariant manifolds that attract every trajectory at an exponential rate. 
They play an important role in the study of the long-time behavior of solu- 
tions, since through them the dynamics of a large system can be described 
by a finite dimensional system. To be more specific, the dimension of the 
state space is reduced by projecting the system onto the inertial manifold 
once the existence of the manifold has been proved. For certain dissipative 
nonlinear systems, the inertial manifold can be shown to exist and it expo- 
nentially attracts solution orbits [16] . In the case where the existence of the 
inertial manifold is unknown, approximate inertial manifolds are introduced 
[TBI [Ml [35l [29] . Several authors considered the construction of approximate 
inertial manifolds |35j as well as numerical simulations using these manifolds 




Inertial manifolds have also been considered for stochastic systems |6l 
[71 [Hi [131 [H]- III one such study. Da Prato and Debussche [9] introduced 
the concept of a stochastic inertial manifold for an abstract stochastic evo- 
lutionary equation in a Hilbert space H (with scalar product (•,•) and the 
induced distance d{-,-)) 



where A is a linear operator, / and g are two nonlinear functions, and Wt 
is a Wiener process taking values in another Hilbert space U. 

Although a theoretical framework for stochastic inertial manifolds has 
been set up [G] [TH [H], it is desirable to efficiently approximate stochas- 
tic inertial manifolds and perform simulations on them. Recently, Roberts 
introduced a normal form transformation of stochastic differential systems 
when the dynamics contains both slow modes and quickly decaying modes 
|33j . in which algebraic techniques were used. This is intended for reduction 
of finite dimensional systems and it might be used to test the numerical 
techniques that are proposed for the infinite dimensional case and also to 
inspire future development. 

In this paper, we introduce a numerical scheme for simulating the stochas- 
tic inertial manifold of a stochastic evolutionary system with multiplicative 
noise, which includes stochastic differential equations (sdes) and stochastic 
partial differential equations (spdes). By projecting to a countable basis. 




[I71[2ll[2l]. 



{ 



du + {Au + f{u))dt = g{u)dWt 
u(0) = no. 
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an SPDE could be converted to an infinite dimensional system of SDEs. The 
main idea is to solve a coupled backward-forward system of sdes where the 
backward part is finite dimensional, but the forward part is either infinite 
or high dimensional. The forward Picard type iteration scheme [H [5l [12] 
is performed on the backward part, and an Euler discretization scheme is 
applied to the forward part. The graph for the stochastic inertial manifold 
is consequently obtained. Two examples, one for a system of sdes and one 
for an SPDE, are presented to illustrate our backward-forward approach. 

The rest of this paper is organized as follows. Section 2 formulates the 
problem and briefly reviews the analytical results [9]. Section 3 introduces 
the numerical scheme and performs the error analysis. An approximation 
procedure is discussed in Section 4. Finally two numerical examples are 
presented in Section 5. 

2 Problem formulation 

We consider a stochastic evolutionary system in a Hilbert space H of the 
form 

du + Audt = f{u)dt + g{u)dWt, (2) 
where the following conditions hold: 

1. Linear part. 

j4 is a self-adjoint operator in H with eigenvalues 

Ai < A2 < • • • < Afc < • • • < Aj +00. 

2. Nonlinear part. 

f : L>(^°) ^ H and g : D{A"') ^ L^(L>(A^)), for some a,/3 G 
[0,1), are globally Lipschitz and bounded, with Lipschitz constants 
Lf and Lg respectively, i.e., for any u,v £ D{A°'), 

I f{u) - f{v) \h< Lf \u-v \d{A'^), 

\ g{u) - g{v) \h< Lg\ u- V \d(^A'-^ . (3) 

When / and g are only locally Lipschitz, but the corresponding deter- 
ministic system has a bounded absorbing set in an appropriate state 
space, it is possible to cut-off / and g to zero outside a ball containing 
the absorbing set. In this case the modified ("cut-ofF') system has 
globally Lipschitz drift and noise intensity. 

3. Noise part. 

Although the Wiener process Wt may take values in a Hilbert space U, 
to be specific, here we just consider Wt to be a two-sided one dimen- 
sional Wiener process, defined in a probability space (fi, J^, P), and 
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adapted to a filtration € M; note that Tt is a two-parameter fil- 
tration [2] or a one-parameter filtration starting at time — oo instead 
of at time [11]. More precisely, for each t E M, 

:= ys<tH. (4) 
jr* := a{W{u),s <u<t), 

which is the information generated by the Wiener process W on the 
interval (— oo,t]. We denote J-^^o -^t simplicity here and hence- 
forth. 

The two-sided Wiener process Wt is defined in terms of two indepen- 
dent Wiener processes Wt and Wt (i > 0), as follows, 



Wt 



I Wt, if t > 0, 
[W^t, ift<0. 

The adaptedness means Wt is measurable with respect to Ft for each t. 

Remark 1. The two-parameter filtration defined in @j requiring C J-^', 
for s' < s < t < t' is consistent with the well-known filtration for posi- 
tive time. Indeed, {J-t}t>o is {^}t>o in the two-parameter setting. The 
only difference between one-parameter and two-parameter filtrations in the 
above setting is their starting time. Since the filtration specifies how the 
information is revealed in time, the property that a filtration is increasing 
corresponds to the fact the information is not forgotten. However, the gener- 
alization of filtration from one-parameter to two-parameter, while maintain- 
ing the property that a filtration is increasing, results in technical difficulties 
in constructing invariant manifolds since a backward SDE is encountered. 
Overcoming this difficulty will be discussed in detail in Remark\^ 

As discussed in O [I3l [HI [15] , it is appropriate and convenient to con- 
sider the canonical sample space, by identifying sample paths of the Wiener 
process Wt with continuous curves (passing through the origin at t = since 
Wq = 0). Namely, a sample path is now a point in the space C(M,]R) of 
continuous functions: Wt{uj) = uj{t). Therefore the sample space is taken to 
be 

n = {uj € C{R;R) : w(0) = 0} 

and P is taken to be the Wiener measure. This is analogous to the situation 
of dice-tossing, where we take six face values, 1, 2, 3, 4, 5, and 6, as samples in 
the canonical sample space = {1,2,3,4,5,6}. When we "toss" a Wiener 
process Wt, we see continuous (but nowhere differentiable) curves as "face 
values" or samples. 
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The Wiener shift 9t is defined as a mapping in the canonical sample 
space 0, for each fixed t G M, 

uj ^ Cj such that 9tuj{s) = ub{s) = uj{t + s) — io{t), s, t G M. (5) 

Remark 2. The Wiener shift defined in ^ is a measure preserving trans- 
formation, i.e. 

F{A) = F{9^^A)yAe 
where P is the Wiener measure and ¥{A) = ¥{9tA) is implied. 

By a simple calculation, we see that 9o = Id (the identity mapping in 
O) and 9s+t = 9s o 9t. Hence the Wiener shift is a deterministic dynamical 
system (or a flow) in il.. The above equation ([5]) means that 

Ws{9tu;) = Wt+s{uj) - Wtiu) « dWtioj). (6) 

Thus 9t is closely related to the noise in the stochastic system 1^ and is 
often called the driving flow. The solution mapping satisfies the property 

m 

S{t, s; oj)x = S{t — s, 0; 9sio)x, x £ H, t > s, F a.s. (7) 

Definition 1. A stochastic inertial manifold for ^ is a random family of 
manifolds {M{uj)}uien which is measurable with respect to Fq and satisfies 
the following three properties: 

1. Each realization of is a deterministic manifold: M{ui) is a Lip- 
schitz (or smooth) manifold for F— almost all lo £ i}. 

2. Invariance: 

S{t,0,uj)M{u}) = M{9tuj), P a.s. 

3. Exponential attraction: for any x £ H 

limd{S{t,0;oj)x,A4{9tCjj)) = 0, exponentially in L'^{Q), 

where u{t,s;u}) is the unique solution for ([2]) defined for t £ [s,+oo) such 
that 

u{s, s; uj) = Us 
and S{t, s; u) is the solution mapping 

u{t,s;uj) = S{t, s;uj)us{uj). 
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Remark 3. 1. M{lo) is a "random variable" mapping all the samples 
(continuous functions) in the canonical sample space into the space 
of deterministic Lipschitz (or smooth) manifolds [1^, and thus is mea- 
surable with respect to J-q. FigureUl explains the meaning of random 
family of manifolds {M{uj)}^(zq heuristically. 

2. The stochastic invariant property can be considered as a natural gen- 
eralization of invariance property in the deterministic setting. In the 
deterministic case, S{t,0)^A = ^A, t G ]R+, which can be seen as 
S{t,0)A4{uj) = since there is one and only one sample in such 

setting. In the stochastic setting, the invariance is in the sense of prob- 
ability. More precisely, under the system evolution or solution map- 
ping S, if we start from somewhere on the manifold A4{u), then after 
time t, we will stand on manifold M{6tU)) whose likelihood of occur- 
rence is the same as A4{uj). This is implied by the measure-preserving 
property of Wiener shift 6t, i.e. ¥{{oj}) = ¥{{6tOj}), see Remark\^ 
From now on, we do not distinguish "uj" and "6tU}" when we say "a 
fixed sample". 



Figure 1: Sketch of a stochastic martial manifold: Three random samples (uj, uj,Cj) 
on the left, and the corresponding three realizations (Ai{Ld),Ai{Ld),Ai{uj)) of the 
manifold on the right. 



Recall that the construction of an inertial manifold in the deterministic 
case amounts to finding a graph above an eigenspace of the linear operator 
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A. By analogy, we take a projection P to a finite dimensional eigenspace 
H~^, and look for a function $ from H'^ = PH to H~ = [I — P)H whose 
graph is invariant under the evolution of the stochastic system ([2]). The 
analytical foundation of our numerical method is the combination of two well 
known methods for constructing deterministic inertial manifolds. One is the 
Lyapunov-Perron method, and the other one is the graph transform method 
[8]. Roughly speaking, the Lyapunov-Perron method looks for solutions 
of the original equation whose components on {I — P)H are bounded for 
negative time. The graph transform method is to let an inertial manifold 
(typically, the flat manifold PH) evolve under the system evolution and to 
verify that the image of PH at time t is a graph which will converge to an 
invariant manifold as t — )• oo. 

As introduced by Da Prato and Debussche [9] , we reformulate ([2]) into a 
backward part and a forward part for time t G [— T, 0]: 

u:=X + Y, 

and 

' dX + AXdt = Pf{X + Y)dt + Pg{X + Y)dWt, 

dY + AY dt = Qf{X + Y)dt + Qg{X + Y)dWt, ^ ' 

. y{-T) = 0, 

where P is a projection from H to the eigenspace spanned by the first k 
eigenvalues Ai, • • • , of yl, and Q := I — P. Here k is determined by the 
eigenvalues and Lipschitz constants for the existence of stochastic inertial 
manifolds [U [9]. Hence X is of finite dimension k and Y is of infinite 
dimension. 

The problem ([8]) involves a backward stochastic evolutionary equation. 

Remark 4. At first glance, we might want to solve the equations ([8]) whose 
unknown is the pair {X,Y), in the interval [— T, 0]. However, this type of 
problem does not have solutions in general. For the existence and unique- 
ness of solution of SDEs, in addition to the usual requirement as in the case 
of ODES, it is also necessary for the solution to be adapted to the filtration 
generated by the noise. Since the filtration is a collection of fields, F := 
{. . . , T-T , J^-t, J^o, J^t, J^T , ■ ■ ■}, J^t C Tt+i, if we use the usual backward in- 
tegration method, i.e., finding the solution at time t by making use of the 
solution at time t-\-l, then this t— solution is Tt+i — measurable but not nec- 
essarily Tt— measurable, which violates the definition of solution for SDEs. 

In order to overcome this difficulty, the terminal value problem of SDE is 
reformulated in such a way as to allow a solution that is {J-t}t<o— adapted. 

By Proposition 3.1 introduced by Da Prato and Debussche [9], for every 
Xq that is J-Q— measurable and square integrable, there exists a unique triple 
{X,Y,Mt) such that 
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1. X : [— T, 0] I— )■ PH is mean-square continuous and adapted, 



2. Y : [—T,0] I— )• QD{A°') is mean-square continuous and adapted, 



3. Mt is a square integrable martingale with values in PH, 

and {X,Y,Mt) solves the following combined backward- forward stochastic 
system, for time t E [—T, 0] 



dX + AXdt = Pf{X + Y)dt + dMt, 
X{0) = Xo, 

dY + AYdt = Qf{X + Y)dt + Qg{X + Y)dWt, 
Y{-T) = 0. 



(9) 



where 



Mt 



Xt-Xo- / AXsds + / Pf{Xs + Ys)ds. 
Jt Jt 



Note that the solution of this system Q sits on the interval [—T, 0], the 
first k components of u, i.e. X, travel backward from to —T, and the 
remaining infinitely many components of u, i.e. Y , travel forward from —T 
to 0. 

Figure [2] is a schematic illustration of this backward-forward method. 



[lo] X QH 




Figure 2: Schematic illustration of the backward- forward approach for approxi- 
mating a stochastic incrtial manifold. 



For each fixed T > 0, define a mapping via 

y(0,w) := <^t{Xo{uj),uj), uen. 



(10) 
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Da Prato and Debussche [9] showed that the hmit of {^t)t>o as T — t- oo in 
L^(ri) is the function <I> whose graph above PH is the inertial manifold 

M{uj) = {u = X + Y : Y = ^{X,uj),X € PH}. 

In the next section we discuss how to numerically simulate ([8]) and then in 
Section 4 we approximate the mapping $ and thus obtain an approximate 
inertial manifold A4. 

3 A numerical scheme 

In this section, we devise a numerical scheme to compute the solution of 
the backward- forward stochastic system ([8|). Our scheme is inspired by [1]. 
Some related references are [5l |T2] . 

3.1 Main idea 

Literally speaking, the backward-forward numerical iteration scheme for Q 
works as follows: We start the first iteration by setting Y to be zero (flat 
manifold) for the entire time interval [— T, 0], and X to be zero on [— T, 0) to- 
gether with the terminal X value to be any J^o-™easurable random variable. 
We obtain the X trajectory backward in time, by using future information 
of (X, Y) at the previous iteration, and then we generate the Y trajectory 
future in time using past in time information of {X, Y) at the previous iter- 
ation. The iteration is stopped when the distance of two consecutive (X, Y) 
trajectories is less than a preset tolerance; we then use the terminal Y value 
at that iteration to approximate one point on the manifold We will 
illustrate the method in more detail in Section [H 

Before getting to the backward-forward approach, we need the following 
preparatory work. 

Let {ei,e2, . . .} be an orthonormal basis for H. In numerical analysis, 
we approximate the Hilbert space (could be infinite dimensional) by a (/c + 
i)— dimensional subspace, and project ([2]) into this subspace, i.e. 

du^+' + Ak+iu''+' dt = fk+i{u^+')dt + gk+i{u^+')dW. (11) 

In theory, / could be infinite or a big natural number. The above projection 
is usually done by the Galerkin method. We denote the numerical approxi- 
mations of X, Y by X, y, respectively. Then PH = R*^ and QH = M', 

= x + y, X eR'' and y G M', 

and 

' dx + Sxdt = ^'^^^{ei,f{X + Y))dt + dMt, 
x{0) = xo, 

dy + Uydt = E^=^+l{e^, f{X + Y))dt + ES+i(e., + Y))dWt, 
y y{-T) = 0, 

(12) 
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where Sx := Ei=l{e^, Ak+iu''+') ,Uy := ES+i(ei, 
For each fixed T > 0, define a mapping via 

y{0,uj) := ^t{xo{uj),uj), u € n. (13) 

The hmit of {^t)t>o as T — >• oo in L'^{Q,) is the function and its graph 
above is the (approximate) inertial manifold M 

M{uj) = {u = X + y : y = ^{x,uj), x G M*"'}. 

Therefore, our goal of constructing the stochastic inertial manifold M is to 
numerically solve (jl2p in order to obtain for a large T > and for a 
number of sample w's. 

For convenience, we denote 

k 

Y,{e^J{X + Y)) ■.= h{x,y), 

i=l 
k+l 

{e,J{X + Y)) ■.= f2{x,y), 

i=k+l 
k+l 

{ei,g{X + Y)) ■.= g2{x,y) 

i=k+l 

for the rest of this paper. In the following, we also denote x{t) by xt and 
y{t) by yt. 

In principle, the solution of (|12p can be obtained as the limit of a Picard 
type iteration (^("^y^"')) as introduced by Da Prato and Debussche [9]. To 
be more precise, {xf'\yl^^) = (0,0), and (x|"\y|"^) is the solution of the 
following iteration scheme 



,{xl''-'>,yr")ds 



+ j% e-^^'+^~^^g2{xf-'\yt-'^)dWs, 

where t G [— T, 0]. The conditional expectation given Tt is introduced to 
guarantee that the solution is adapted, i.e., measurable with respect to the 
filtration Tt , as is done in the theory of backward SDEs |3H [32] . 

Our goal is to find ^>t via looking for y(0,a;) for given xq since y(0) := 
^T{xo-,i^)- We now introduce a time discretization of the above iteration. 
Note that for the backward part x, the conditional expectation is still in- 
volved; for the forward part y, we use the Euler-Maruyama scheme in the 
discretization. In the numerical computation, every simulation corresponds 
to one sample, however, we cannot specify which sample we have actually 
chosen, but we do know that it is a validated sample in the sample space Jl. 
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Suppose h := T/N, ti = —T + ih, i = 0,1, . . . , N , taking Wt to be 
one dimensional Wiener process, and denoting AWi := Wj+i — Wi, then a 
time discretization of (|14p is 



h 



In Section [U we will explain how to implement this scheme. Note that 
we could also use the available {x[^\y^J^), instead of {x[^ ^^^ytj ^^)' at 

the current iteration (for j = 1, • • • , i — 1) to calculate in the above 
scheme (fTSll . 

We now show the convergence of the above time discretized Picard iter- 
ation scheme (fTSl) . 

3.2 Convergence of the numerical scheme 

We now prove the convergence of the numerical scheme (jlSp devised in 
the last subsection. Namely, we prove the convergence of {x["'\yl-'^^) to 
{xt,yt) in a certain sense (see Theorem [1] below). To this end, we estimate 
Xj"^ — x[°°\yl^^ — y^'^ and — xt,y^'^ — yt in Lemma [Hand Lemma O 
respectively. 
Define 



(oo) 

(°°) n 

Vto = 



I = e-^'[ytl + IL e^^^-'-^^hixtlytl)ds + e^i^-^^-^^ g^ixiZUtbdWA- 

Recall 

=0, < f < iV - 1, xtj^ =xo 

yif =0, 0<j<iv 

xt'^ = iE{(xo + Ef=^^Nf - hi4:\yi:^)]h) I 

[ yf'^ = e-^'iyt''^ + It^ e^(^-'^-^)h(xi^^,yl^^)ds + f^^ e^(^~^^-^) 9,{xP_^,yi:l)dWs]. 

Lemma 1 (Iteration error). Let the Lipschitz condition ([3|) be satisfied. 
Assume that 

I /i(0,0) I + I /2(0,0) I + I 52(0,0) I + I xo |< i^. 
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and 



M- 



1 := max{2L^/i, 6e-'^^^{Ljh + L^)} ^^^^ y^M{T < 1. 



Then the following inequality holds: For all < i < N, 



max 

0<i<N 



E I X 



in) 



X 



(oo) 1 2 



+IE I yl 



in) 



(oo) 1 2 



<c2(-+c7i/ir, 

where Ci depends on Lf,Lg,U,T and C2 depends on Lj, Lg, K,U,T 
Proof. First note that 



X 



(n+l) 
M) 



E{x 



(n+l) 



\Sx 



in) 



We now estimate 

< e(e| I - x'^^ + [sx^^ - &sr'^ - (/i(4r^2/ir^) - hi4r'\yt'^))]h\'\ ^u]) 

(Jensen's inequality) 

= e{ I - <\ + Nr^ - ^4""'^ - (/i(4r^2/ir^) - p } 

= e{ I - <\ I' + 1 - sxf^ - {fi{4:\yP) - fMr'\yt'^))] h \' 

+ 2| (xf+^)-x(^J^/^|| [Sxi:^-Sx[:-'^-if,ix^:\y[f)-fMr'\yt'^))]Vh 

< e{ 1 - <\ p + 1 Nr^ - ^xSr'^ - ) - fMr'^^yf^mi' } 

+ 27e[| xS:;t'^ - <\ p I Nr^ - ^4""'^ - ) - /i(xir^\yjr^)))]^/. i 

< Ie{ I 4Z'^ - <\ P + I [SxSf) - Sxf^ - (/i(xf ), yjf)) - /i(xlr '\ yf-'^))] h \' } 
+ r/.E I x'C^^ - xfi P +T-'hE I 5xf ) - 5xlr'^ - mx^\yi:^) - Mxif,yt'^)) p 
(Young's inequality ab < — H — — where T > will be choosen later) 

<(i + r/.)(E|xSf+^)-x;ri p 

I h + r~^ 2y I An) J"-l)|2, /^+r2^.2|p| M (n-1) ,2 , A 
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Since x 



tN 



xi = xq, by iterating the last inequality, we obtain 



E[| X 



(n+l) Jn) |2i 



N-1 



<L}ih + T-')[Y,^ 



(n) (™-l) 1 2 

yl - yl I 



Af-l 



(n) |2 



0=1 



h]. 

(16) 



For the forward SDE part, since in the numerical scheme y is only of finite 
dimension, we consider y to be one-dimensional for simplicity, and estimate 



E \e 



-Uh 



U-i 



+ I e 



-2Uh 



-2(7 /i 



-2Uh 



IE 1 yt;"' 


(n) 


2+E| r e^(^-*-^)(/2(xt\ 


,2/t\)-/2(-t";'^?/£;'^))^« 


+ E| / 


^U{s-t 

1 






E 1 .tr^ 


(n) 




1)^ 


+ hELli\ 


(n) 

— 






IE 1 yt'^ 


(n) 


2+(2/,2^2^2/iL2)E|xt\- 


^.("-l) |2 



<6e-2^'^(/iL^ + L2)( ^E 



X 



(n) |2 



i-1 



/i + ^E 



(n) |2 



%\ -yi 



(17) 



where the last inequality is by iteration. 

Combining ([16]) and ([IT]), and letting Mi := max{L^(/i+r-^), 6e"2^''(L^/i+ 
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Lg)}, we have 



max 

0<i<N 



I (n+1) (n) 1 2 



< MiT max 

0<j<Af 



E I x^"'^ — Xi 



+E I yl 

|2 



|2 



, IT, I (n) 1 2 



< (MiT)" max 

0<j<Af 

= (MiTV max 

0<i<Af 



E I X 



(1) 



,(0) |2 



+E I yl'^ - yi: 



,(0) |2 



E I xll> |2 +E I 



Define (MiT)" maxo<i< 



j<7V 



E I p +E I yil^ |2 



by On- Then 



the series a„ converges when y A/iT < 1. Therefore, for all < i < A^, 

{{x[^\yt^^)}n is a Cauchy sequence and thus converges to (x 
the mean square sense. To be more precise, 



(oo) (oo) 



m 



max 

0<i<N 



E I X 



(oo) 



X 



(n) 1 2 



+IE I yl 



(^) 



(n) 1 2 



■ u=n 



E 

(MiT)"maxo<i<7v 



E I xSj) |2 +E I yj^'^ 1 2 



(1- VAfir)2 



Choosing T := 1/h, Mi = Mi as in the assumption, thus \/MiT < 1. 
By taking MiT := ^ + Cih where Ci depends on Lj, Lg, U, T and Cih < i, 

since maxo<i<ArE | x^^'^ p +E | y^^^ p< C3 as shown similarly as Lemma 8 
in Bender and Denk [3] and C3 only depends on Lf,Lg,K, we have 



max 

0<i<N 



E I x 



(~) 



{") 1 2 



+IE I yi. 



yi 



{") 1 2 



<c2(-+ci/ir, 



where C2 depends on Lj, L^, K, U,T. This proves the lemma. 



□ 



In the following, C will denote a generic positive constant, independent 
of i and n, that may take different values from line to line. 

Lemma 2 (Discretization error). Assume all the conditions as in LemmaU\ 

are satisfied, then 



sup 

-T<t<0 



E\xt-xt^ \^ +E\yt-yt^ \^ ] < Ch. 
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Proof. Note that 

For ti^i < t < ti, we have 

4zl = nxt^ I ^u-A + [sxtl - fMT.UrJM - ^-i)- 

By the Martingale Representation Theorem 



= nx^r^ 



we have 



(oo) (oo) 
— . 

which yields 



By Ito's formula, we have 



Let 



6xt := xt - x[°°\ 6Zt := Zt - Z^°°^ 



t ' 



and 5ft := (Sxt - fi{xt,yt)) - (^H - fiixtlytl))^ 
then define 

:= E I fet P + ^ ' E I SZs P - E I fet. ^ 
^ ' E[25x,5/,]ds 
<E[C^*'|fe, 1(1 + 
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Since 

E\xs- xt,_j p = E 
< 2 



r {-SXr + fl{Xr,yr))dr+r ZrdWr? 
[E \ / {-SXr + fl{Xr,yr))dr\'^ +E\ / ZrdWrl"^] 

f'S f'S 

<2[{s-ti_i)E {-Sxr + fiixr,yr)fdr + E Z^dr] 

<2[{s-ti.i) r [E(-5x, + /i(x,,y,)-/i(0,0))2+E(/i(0,0))2]dr 
+ E / Z^dr] 

< 2[(s - ti_i) I E(C7x2 + 2Lf{xl + y^) + C)dr + E f Z^dr] 

< Ch, 

then we have 

E I X, - p < 2[E I xs - xu_, |2 +E I 5xu_, p] 
<C{h + E\ 5xt^_^ |2) 

for ti-i < t < s < ti. On the other hand, 

E\ys- yu_, \^<C + CE\ f e^(^+^)/2(x„ y,)(ir + f e^^'+^^ g2{xr, yr)dWr P 

<C + C[{s-ti.,)E r e2^(^+^)/|(x„y,)dr + E r e'^(-+''^ gl{xr,yr)dr] 

< Ch, 

and note also that y^°°'^ agrees with y at each grid point, i.e. yl'^l = y^^i, 
thus 

Ely.- ytl I' < 2[E | y, - yu_, p +E | p] 
< C(/i + E I (5yt,_, |2) 
= C/i. 

By the definition of At, 

At< aE \ 6xs l"^ ds + — [/i + E | Rds 
Jt a Jt 

< [\e\6xs ds + -[h^ + hE\ 5xt^_, 
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and 

E I fet P < E I fet P + / El 6Zs 1^ ds 



< J^\e\5xs ds+(E\ 6xt, 1^ +^[^^ + hE \ (5xt,_, 1^]^ . 

(18) 

Let i?i := E I ^ +^[/i^ + /iE | 5xt^_-^ p], then by Grownwall's inequality, 

E I dxt |^< i?ie^«'"'^* < Bic"'^. Plugging it in the second inequality of ([T8]), 
we have 

E\6xt\^ + J^'k\6Zs f ds < Bi{l + ae''''h). (19) 
By taking t = tj-i, we have 

E I 5xt_. 1^ + / E\5Zs 1^ ds < (1 + C/i)(E | fe*. ^ + hE \ 5xt_, p 

Ju-i " 

and if we choose a ^ C, 

E I p + / ' E I p < (1 + (:7/i)(E | 6xt^ p +/i2). (20) 

Iterating the last inequality and recall that E | Sx^^ p= 0, we have for 
sufficiently small h, 

E I Sxt^_, |2< (1 + Ch)^{E I 6xt^ |2 +h) < Ch, 

which yields Bi < Ch, and by (I19p . we get for all t G 

E I 5xt p< Ch, 

and the right hand side of the inequality does not depend on t. Therefore, 
sup E I 6xt P= sup E I - xt P< Ch. (21) 

-T<t<0 -T<t<0 
For the forward part y, we have 

E I yt^ -yt\' = e( f e-^(*-)[/2(x(r;,yf J) - h{xs.ys)]ds 

t , , , \ 2 



< 2/iE / e-2^(*-)[/2(4°°;,yJ~;) - /2(x„ 2;,)]2ds 

+ 2E f e-2^(*-^)[,2(xtl.?/H)-52(x.,y.)]2cis 

< [C(/i + E I 6xt,^, |2) +C{h + E\ 6xt^_, \^)]{2h^C + 2/iC) 

< C/^^ (22) 
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where the last inequahty is due to (pO|) . Consider (f2T]) and (p2]) . we have 
sup E I - |2 +E I - yt \^< Ch. 

-T<t<0 

This proves the lemma. □ 

Now we are ready to state the convergence theorem. For t G 
define xj""^ = x[^l^ and yj"^ = y^^^"^^^ . We can also define Xj"^ and as the 

linear interpolation among Xj^^'s and among y^^^^^'s, respectively, and the 
following result also holds. 

Theorem 1 (Convergence). Assume that all the conditions in LemmaUl are 
satisfied. Then 

sup (E\xt- xl") |2 +E\yt- A< C{h + + C/i)"), 



-T<t<0 



-2 



where C > wiZZ denote a generic constant depending only on the data 
Lf,Lg,K,U,T. This implies the convergence of the numerical scheme ([T5]) . 
as /i — 7- and n — )■ oo. 

Proof. Note that 

E|xt-x|"^ |2+E|yt-2/f) I' 
<E I xt - xt^_, I +E I xt^_, - x\J^ I +E I xl_[ - xj,_\ I 
+E I y, - y,^_, |2 +E | y,^_, - y\Z\ ? +IE | vtl " V^l l' 

<c{h + {\ + chT) 

by regularity of the true solution as well as Lemma [J and Lemma [2j Taking 
supremum on both sides of the above inequality, we have 

sup ("e I xt - xt^ |2 +¥.\yt- yj"^ \A < C{h + + ChT). 



-T<t<0 



□ 



4 Approximation of the stochastic inertial mani- 
fold 

Now we approximate the graph $ for the inertial manifold M. Recall that 
$ is approximated via 

y{0,uj) := ^>r(xo,w), u e (23) 
18 



When T is sufficiently big, ^ ^- So we need to evaluate 1/(0, a; ) (which is 
y(0) in Section 3). To this end, we need to compute, step by step, x[^\y[^^ 
in (jlSp . For the backward part x["'\ the conditional expectations E[- | Tt..] £ 
L?'{J^ti) in (|15p will be approximated by their orthogonal projections Pi on 
finite dimensional subspaces Aj of L^(J"f.), where L'^{Fti) contains all the 
functions in Lp'iO) that are J^j. -adapted. 

in) 

Indeed, instead of computing follows 



xt^ = E< 



|L-^(sx(;-^) + /i(. 

3=i 



,(n-l) ^ (n-1) 



we will compute its orthogonal projection x\. : 

N-l 



X 



in) 



j=i 



(n) 

(n-l) ^(n-1) 



. yh 



))A, 



have 



Denoting a basis of the projection space Aj by {r]\, . . . we then 



X 



D{i) 

(n) _ (n) 



(24) 



Set 7?* := (r/^, . . . , ??})(j))' and a-"^ := (a-"^'', . . . , "-^(i))', with prime denoting 
matrix transpose, the calculation for the backward part x^. boils down to 
two aspects: one is the basis ry\ the other is the coefficient . 



Remark 5 (Basis). As in j^Tf , a complete orthonormal basis of L {U,J^t,^ 
is given by the Wick polynomials {T*(^);a G J'}, 

where Hai{-) are Hermite polynomials. Here 

oo 

J ■= [a = (aj, i > 1) I Oj G {0, 1, 2, . . .}, I a 1= ^ Oj < oo}. 



(t) := / mi{s)dWs, 
Jo 



and mi{s),i = 1,2,..., are a set of complete orthonormal basis in the Hilbert 
space L^([0, t]). One choice of the basis elements are the Hermite polynomi- 
als of Brownian motion Wt ■ In fact, 
1 



= (1,0,0 

a2 



(2,0,0 



m) = H^iumoiumom)) ■■■ = Hiim), 



(n, 0, 0, • • • ), r*„(^(t)) = H2iUt))HoiUt))Ho{Ut)) ■■■ = ^2(6(t)). 
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Since ^i(i) = jl^mi{s)dWs = j^ldWa = Wt, the basis elements "rf^', 
HfiiyVt^) which are adopted in our numerical implementation. 

In principle, the coefficients af''^ are calculated as follows: 

N~l 



s are 



OL. 



in) 



j=i 



(25) 



where /3j := I E [77^77*] J are the inner- product matrices associated 

V J p,q=l,...,D{i) 

with the basis. 

In practice, /3'-s are computed by their simulation-based estimators, such 
as Monte Carlo least squares estimators used here [3]. To this end, we 
are assuming to have r (sufficiently large) independent copies {xWi,\r]li), 
A = 1, . . . , r, of {Wi,r]^^). Here the index A denotes copies. Then 



ft = ( nvi%. 



( Eri\r]\ Er]\r]\ 



p,g=l,...,D(i) 



Er 



\ 



is replaced by its Monte Carlo simulation /3j 



ft 



1 



and is further rewritten as 



ft 



{Ai)'Ai = -f X^A?/aA?/b), a,b 
^ ^ A=l ^ 



J2l=i\vlxvl I 



l,...,I)(i), 



where 



= — xrfd 



( 1??2 \ 

2??i 242 ^^(i) 



\ rVl rV2 ' ■ ■ rVD(i) J 

The pseudo-inverse denoted by := {A'A)~^A' is used in the computation 
of ft. 
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The calculation of the coefficients a^"^ is to obtain the backward part 
x[^^ . Note that the calculation of a^"^ also needs the forward part ^\ 
which is calculated through Euler-Maruyama scheme. The following formu- 

(n) 

lae illustrate how to update a\ . 



D{i) many 



(26) 



jA=l,...,r — 



-(n-1) .(n-1) 



(n-l), 



D{i) D{i) D(i) 

fSr^ (n-1) i sr^ {n-l) i 

d=l d=l d=l 



Y^a^'^rvii)', (27) 



(Ay(°))A=i, 



0, 



( \Xq \ 

22;o 



(n-l) 
J 



(28) 

(n-l) .(n-l) 



A 
(29) 



1 



A=l,...,r 



A, 



(30) 



(A^?^)A=l,...,.^(lC^2C^•••-0' 



-;(")V 



+ 52(Axt-/\Ayt;^VAW.] 



(31) 



A=l,...,r 



.(n-l) ^ . ,(n-l) 



are independent copies of x^" cor- 



where axJ" = (i^t 
responding to independent copies of basis functions 
Upon converge, i.e. 



max{E I X 



(n) ^("-1) |2 



(32) 



we have y['^ as our approximation of yo- 

Remark 6. Although we only need the final grid point value ytj^ (recall that 
t]\f = 0) as our approximation of the inertial manifold y^ := <I>r(xo), the 
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intermediate points (xj" , y^" ) are approximated by (xj" , yj" ) as follows 
4:\u^) = Y.^Svk^)' (33) 

d=l 

(34) 

Therefore, this approach of approximating stochastic inertial manifold also 
provides a way of solving backward-forward stochastic differential equations. 
As in equation (|33p . x[^\uj) 's are approximated by its orthogonal projection 
x^^^\uj) := Yld=i '^td ''Idi^) ' where the randomness comes from the basis 

(n) 

functions rf^ rather than the coefficients ^ . The fact that a 's are deter- 
ministic can be seen from (j25p . since a's are expectations. In the numerical 
simulation of a's, we utilize copies of{x,y)'s (different copies corresponding 
to different sample paths oo) to calculate the expectations. 

All the conditions required by Bender and Denk in [1] are satisfied in 
our case, so the error analysis results for the Monte Carlo simulation also 
apply here. 

Figure [3] demonstrates the procedure of computing stochastic inertial 
manifold. 



Stochastic systems (1) 



Galerkin projection (10) 



Backward-forward stochastic system (11) 



For each Xq, solve (11) to find yo 



Inertial manifold is plotted by linear interpolation 



Figure 3: Proccdmv for computing a stochastic inertial manifold. 



Specifically, the computation is achieved in the following way as shown 
in Figure U) We begin with the flat manifold, i.e. let all the y-values be 
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zero. In principle, we could take any other acceptable initial manifold and 
let it flow forward, and at the limit it would also approach the desired man- 
ifold. We thus start our iteration by letting the initial guess y be zero at 
each tj, see (|28p . In order to avoid nestings of conditional expectations, we 
also let initial guess x values to be zero except at terminal time tjv, which 
corresponds to letting all but the final set of coefficients a to be 0, see (j26p . 

(n) 

Recall that we approximate each xj, value by its finite dimensional orthog- 
onal projection as in (|24p . so indeed we want to calculate a^"^- The terminal 
value of X is set to be a J-q— measurable random variable. At each iteration, 
in order to update a."^-^ (pO|) . we use copies of • • • ,a a^j" as well as 

^yti 1^^' " " " 1^ ('^ = 1) ■ ■ ■ 5 ^) from previous iteration, and then gener- 

ate copies of x[^l^ by virtue of copies of basis functions r/'s, see Copies 

of Ut^^-^ are reproduced by an Euler-Maruyama scheme as in (j3ip . in which 
different copies correspond to the sample w's that are already chosen in basis 
functions 77's. This procedure is repeated until {x, y) converges in the mean 
square sense ()32p . We finally acquire the terminal value of y at the stopped 
iteration as the approximation of yo- 



y 



Stochastic inertial manifold ya = if>{x,Lj) 





;r(2){-T) x-<i'(-T) xa{ijj) 



X 



Figure 4: Schematic diagram of the numerical iteration. 
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5 Examples 



In this section, we test our backward-forward numerical scheme in two ex- 
amples. One is a system of SDEs, and the other one is an spde (which is 
converted to a system of sdes). 

Example 1: A system of stochastic ordinary differential equations 

dXt = {-aXt - XtYt)dt, (35) 
dYt = {-Ytil + 2Yt)+ + X^)dt + aYt o dWt, (36) 

where Wt is a scalar Wiener process, a and a are real parameters, o indicates 
the Stratonovich interpretation of the noise term and = max(/, 0). 

Figure [S] is the phase portrait for the deterministic counterpart of Ex- 
ample 1 ((T = 0). 
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Figure 5: Slow manifold for the deterministic counterpart of Example 1: cr = 
and a = 0.1. 

Take o" = 0.1 and a = 0.1. Figure [6] shows several sample solution paths 
of the SDK system (f35D - p6D . 

Roberts [33j introduced a normal form transform method for stochastic 
differential systems with both slow modes and quickly decaying modes. The 
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X 

Figure 6: A few solution paths of Example 1: a = 0.1 and a — 0.1. 

(approximate) formula for the slow manifold (a type of inertial manifold) of 
Example 1 obtained via his method 0. is 

y^x^ + OAx"^ f e-^^-^^dWr{uj) (37) 

and 

x^x + O.lx I e-^^-^UWr{u:). (38) 

In Figure O we plot the stochastic inertial manifold according to y ~ 
x2 + ^.Ix"^ j_^^e''dWr and compare it with the stochastic inertial mani- 
fold from our backward- forward method for one sample path oj. There is a 
remarkable agreement of the shapes of the stochastic slow manifold obtained 
by the two distinct methods. For four samples in Figure [HI Figure [9] then 
shows the discrepancy between the stochastic inertial manifold realised on 
three realisations oj, Co, Co and that realised on oo via our backward- forward 
approach, respectively, i.e. M{oj)-M{ijj), M{Cj)- M{uj), and M {Cj) - M {uj) 
in red, blue and green color. 

Example 2: A stochastic partial differential equation 

: //www. maths . adelaide . edu.au/anthoiiy.roberts/ sdesm.html| 
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Comparison of SDE slow manifolds constructed via b-f sctieme and normal form transform 
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Figure 7: Stochastic inertial manifold A4{uj) calculated by the backward-forward 
method and by the normal form transform method for Example 1. Online ver- 
sion: Magenta from the backward-forward method and black from the normal form 
transform method. 

We examine how effective the stochastic inertial manifold approach is in 
assisting the simulation of the long-term dynamics of an spde. We consider 
a stochastic partial differential equation 



where v = 0.01, Wt is a scalar Wiener process, and Ito interpretation of the 
noise term is adopted. 

Let us first consider the dimension of the inertial manifold. Note that the 
eigenvalues of the operator vdxx + 1 are —vi^-k"^ -\- 1 with the corresponding 
eigenmodes ei{x) = \/2 sin(i7rx), for z = 1, 2, • • • . Thus the dimension of the 
deterministic unstable eigenspace is 3. In this example, H = L^(0,1). Let 
P be the orthogonal projection to the (deterministic) unstable eigenspace 
i/"*", spanned by eigenmodes ei,i = 1, 2, 3. 

The existence of the stochastic inertial manifold requires the nonlinear 
terms to be the globally Lipschitz. Here, in Example 2, although the drift 
term f{u) := u — v? is only locally Lipschitz, we can prepare the equation by 
replacing / by the cutoff function / which is defined to be / in a bounded 
neighborhood centered at the origin, and otherwise. The details of this 
procedure were described by Da Prato and Debussche [9]. 




(39) 
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First sample uj 



Second sample u 




Third sample uj 




pburth sample uj 



-10 




Figure 8: Four samples uj, oj, Cj and Co for Example 1. 



With the Galerkin projection u = X^j iij(t)ej(a;), the evolutionary equa- 
tions we will be working on become 

dui{t) = v\iui{t)dt + (n — n^, ei)dt + ui{t)dWt, 
p : < du2{t) = u\2U2{t)dt + {u — V?, 62) dt + U2{t)dWt, 
duz{t) = u\zuz{t)dt + (u - u^, ez)dt + U2,{t)dWt, 

q : dui{t) = v\iUi{t)dt + {u — , e/Cjdt + Ui{t)dWt., 

where we only keep four Fourier modes Ui{t)ei,i = 1, 2, 3, 4 for u; U5{t),UQ{t), . 
are simply ignored, as numerical simulations indicate that they are negligible 
in this specific case. 

It is difficult to visualize the stochastic inertial manifold directly. But we 
can plot a "point" p in and the corresponding "point" q = ^{p,uj)e4{x) 
on the inertial manifold, separately. The "point" p may be represented as 
p = X]i=i ^'(Oi '^)^«(^) with coordinates (ui(0,a;), 1*2(0,0;), 1x3(0,0;)). The 
corresponding "point" q on the stochastic inertial manifold is then computed 
through our backward-forward approach. 

We plot three different realizations of a point p versus space variable x 
in the left panel of Figure [T0| and then plot the corresponding point q versus 
space variable x separately in the right panel. 
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■r X 



Figure 9: Discrepancy between different realisations of the stochastic incrtial man- 
ifold in Example 1 via the backward- forward method. Online version: Magenta (top 
left) for the stochastic manifold realised on sample ui, while red (top right), blue 
(bottom left) and green (bottom right) for discrepancies between stochastic man- 
ifold realised on samples Co, Co and vj and stochastic manifold realised on sample 
u. 




Figure 10: Example2-A "point" p — ui(Q,Lo)ei{x)+U2(0,uj)e2{x)+U3{0,uj)e3{x), 
with ui,U2 and U3 being uniformly distributed random numbers, in and the 
corresponding "point" q in the stochastic inertial manifold: Three samples of p 
(left panel) and the corresponding three samples of q (right panel). Online version: 
Corresponding p and q samples have the same color (green, red or black). 
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6 Discussion and conclusion 



In this paper, we have devised a backward-forward numerical approach for 
computing stochastic inertial manifolds for stochastic evolutionary equa- 
tions, including higher dimensional stochastic ordinary differential equations 
and stochastic partial differential equations. This approach is based on the 
stochastic inertial manifold theory of Da Prato and Debussche [9] , which re- 
quires a backward- forward approximation formulation. In fact, our approach 
also provides a stand-alone numerical scheme for solving backward-forward 
SDES. 

Unlike deterministic evolutionary equations, we need to guarantee the 
adaptedness of solution processes with respect to an appropriate filtration 
for stochastic evolutionary equations. This makes the computation in the 
backward part of our numerical approach cumbersome, as it involves sim- 
ulating conditional expectations which further requires a random basis. It 
would be nice to have a numerical scheme that uses only forward simula- 
tions. Clearly, a numerical scheme involving only forward simulations would 
be desirable; such schemes have been devised for deterministic problems (e.g. 
|20^ I36j ) and we are working on adapting them for the stochastic case. 

Acknow^ledgements. We thank Arnaud Debussche, Arnulf Jentzen, 
Edriss Titi, and Jianfeng Zhang for helpful discussions. 
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